Dynamics of tumor-associated macrophages in a quantitative systems pharmacology model of immunotherapy in triple-negative breast cancer

Summary Quantitative systems pharmacology (QSP) modeling is an emerging mechanistic computational approach that couples drug pharmacokinetics/pharmacodynamics and the course of disease progression. It has begun to play important roles in drug development for complex diseases such as cancer, including triple-negative breast cancer (TNBC). The combination of the anti-PD-L1 antibody atezolizumab and nab-paclitaxel has shown clinical activity in advanced TNBC with PD-L1-positive tumor-infiltrating immune cells. As tumor-associated macrophages (TAMs) serve as major contributors to the immuno-suppressive tumor microenvironment, we incorporated the dynamics of TAMs into our previously published QSP model to investigate their impact on cancer treatment. We show that through proper calibration, the model captures the macrophage heterogeneity in the tumor microenvironment while maintaining its predictive power of the trial results at the population level. Despite its high mechanistic complexity, the modularized QSP platform can be readily reproduced, expanded for new species of interest, and applied in clinical trial simulation.


INTRODUCTION
Breast cancer is the most commonly diagnosed cancer with over two million new cases worldwide annually (Ferlay et al., 2021). Triple-negative breast cancer (TNBC), which is defined by its lack of estrogen receptor, progesterone receptor, and human epidermal growth factor receptor-2 expressions, accounts for around 15% of total breast cancer cases. TNBC is considered an aggressive phenotype of breast cancer, and while patients can respond well to chemo-immunotherapy, those with residual disease have high rates of recurrence with approximately one-third developing recurrent disease within three years (Hensing et al., 2020;Schmid et al., 2020). Once metastatic, the disease is incurable and overall survival is substantially lower compared to other types of breast cancer (Lee et al., 2020;Lin et al., 2012). Although chemotherapy remains the mainstay treatment for patients with metastatic TNBC, targeted therapies such as sacituzumab-govitecan have recently been developed (Bardia et al., 2021;Malhotra and Emens, 2020). Biomarker-based treatments are also available, such as for patients with germline BRCA mutation and programmed death-ligand 1 (PD-L1) expression in the tumor (Anders and Carey, 2022). For patients with PD-L1-positive metastatic TNBC, the anti-PD-1 antibody pembrolizumab combined with chemotherapy is a standard firstline treatment based on the results from the KEYNOTE-355 trial (Cortes et al., 2020). For patients with germline BRCA mutations, the poly (ADP-ribose) polymerase (PARP) inhibitors olaparib or talazoparib are recommended (Barchiesi et al., 2021).
Despite the availability of these targeted agents for metastatic TNBC, chemotherapy remains the standard-of-care treatment for PD-L1-negative TNBC without BRCA mutation (Anders and Carey, 2022). Given the relatively low prevalence of PD-L1-positivity and BRCA mutation in TNBC (about 40% and 10-20% respectively), multiple new therapeutic targets are under investigation in preclinical and clinical studies (Pogoda et al., 2020;Torres and Emens, 2022). Currently, the most promising treatment for metastatic TNBC is immunotherapy involving blockade of immune checkpoints, such as PD-(L)1 and lymphocyte-associated gene 3 (LAG3), in combination with other therapeutic agents (Huppert et al., 2022). Immune checkpoint molecules are membrane proteins that act as regulators of the immune system, which suppress antitumor immunity in cancer (He and Xu, 2020). In particular, CD47 is an immune checkpoint molecule that can inhibit phagocytosis of cancer cells by tumor-associated macrophages (TAMs), In the past few years, we have developed and expanded a large-scale QSP platform for the analyses of immune checkpoint inhibitors and bispecific T cell engagers in combination with other agents in non-small cell lung cancer (Jafarnejad et al., 2019;Sové et al., 2020), colorectal cancer (Ma et al., 2020a;2020b) and breast cancer (Wang et al., 2020a. We have also combined the QSP model with a spatial agent-based model of tumor to describe spatial heterogeneity of the tumor microenvironment (Gong et al., 2021;Zhang et al., 2021a). Here, by integrating a macrophage module into our previously published QSP platform (Wang et al., 2020a, we are able to investigate the impact of TAMs on the cancer-immune cell interactions and provide a computational framework to predict clinical response to macrophage-targeted agents based on preclinical data. We aim to demonstrate our calibration and validation steps to integrate the new macrophage module and show that the expanded model retains its predictive power with recalibrated virtual patient population using recently published data. Besides those inherited from our previous study (see STAR methods), new data includes intratumoral cytokine concentration measured directly from TNBC tumor biopsy samples (Autenshlyus et al., 2021), M1/M2 macrophage ratio from omics data analysis (Reynolds et al., 2017), and density of TAMs estimated by digital pathology analysis (Yang et al., 2018).

Integrating the macrophage module into the QSP platform
In our previously published QSP platform for TNBC immunotherapy, modules were built to investigate the dynamics of T cells (i.e., effector T cell, regulatory T cell, helper T cell), antigen-presenting cells (APCs), cancer cells, tumor-specific neoantigens and tumor-associated self-antigens, immune checkpoints, myeloidderived suppressor cells (MDSCs), and therapeutic agents, respectively. Figure 1 shows the interactions between TAMs and the pre-existing model species. As described in our previous study , the model includes four main compartments: central compartment describes circulation of therapeutic agents and immune cells between the bloodstream and the other parts of the body; peripheral compartment represents peripheral organs/tissues with naive T cell maintenance (den Braber et al., 2012); lymph node compartment represents tumor-draining lymph nodes at immediate downstream of a breast tumor, where T cell activation occurs; tumor compartment describes dynamics of cancer cells, activated T cells, ll OPEN ACCESS Figure 1. QSP model diagram The model is comprised of four compartments: central, peripheral, tumor, and tumor-draining lymph node, which together describe cycles of immune activation in lymph nodes, T cell trafficking to the tumor, killing of cancer cells, immune evasion, and antigen release and lymphatic transport . nT, naive T cell; aT, activated T cell; NO, nitric oxide; Arg-I, arginase I; Treg, regulatory T cell; Teff, effector T cell; Th, helper T cell; Mac, macrophage; mAPC, mature antigen-presenting cell. Cytokine degradation and cellular clearance are omitted in the figure.
more rapidly than recruitment, so the entire process is modeled by a CCL2-mediated recruitment of M1like macrophages (Equation S50). Once recruited, TAMs can be polarized toward pro-(M1) or anti-inflammatory (M2) macrophages by T cell-and tumor-derived cytokines. Figure 2A shows the cytokines that participate in the polarization process. IL-12 and IFNg shift TAMs toward M1-like phenotypes while IL-10 and TGF-b shift TAMs toward M2-like phenotypes. The model calibration started by estimating the recruitment and polarization rates of TAMs. The recruitment rate was estimated by the median density of TAMs measured in tumor nests from patients with basal-like breast cancer, a subtype similar to TNBC, by immunohistochemistry using a pan-macrophage marker, CD68 (Yang et al., 2018). The reported 2D density was converted to 3D density using a stereological equation (Mi et al., 2020). Assuming that M1-to-M2 polarization can occur as fast as two days, we estimated the rate of M2-to-M1 polarization to match the M1/M2 macrophage ratio with the median of clinically measured values in basal-like breast cancer. The clinical measurements were extracted from the ISB Cancer Genomics Cloud (ISB-CGC) (Reynolds et al., 2017), where proportions of M1-and M2-like macrophages among immune cells in breast tumors were determined using CIBERSORT (Newman et al., 2015). Notably, the calibration aimed to fit the model-predicted values (i.e., cellular density, M1/M2 ratio, etc.) at the time point when the tumor reached the median pre-treatment lesion size (Table S1), given that majority of the clinical measurements were taken upon diagnosis.
The functions of M1-like macrophages are facilitating APC maturation via IL-12 secretion, and phagocytosis of cancer cells. As IL-12 is the major cytokine for mediating APC maturation and subsequent T cell activation, it is mainly secreted by M1-like macrophages and mature APCs (mAPCs). The secretion rate of IL-12 by macrophages is estimated by the IL-12 level produced by TAMs isolated from solid tissue of human ovarian carcinoma and macrophages during the wound healing process in humans (Clough et al., 2007;Sica et al., 2000), and the secretion rate by mAPCs is roughly ten times higher upon antigen stimulation (Heufler et al., 1996). Rates of cytokine secretion by M2-like macrophages (i.e., IL-10, TGF-b, and VEGF) are estimated similarly (see STAR methods). Further, macrophage-mediated phagocytosis is a slow process that not only takes days to complete but is also inhibited by IL-10 and CD47 and PD-L1, two immune checkpoints on cancer cells (Bian et al., 2016;Gordon et al., 2017;Pan et al., 2020). To model the overall inhibitory effect of checkpoint molecules on phagocytosis, we built a submodule that incorporates interactions that impinge on these checkpoints, including PD-L1 upregulation by IFNg, trans interactions between CD47 and SIRPa, PD-1 and PD-L1, and cis interaction between PD-L1 and CD80 on the surface of cancer cell. We also incorporated checkpoint blockade by specific antibody drugs (e.g., anti-CD47 and anti-PD-L1).
The overall inhibitory effect of PD-1 and SIRPa on phagocytosis is calculated by Hill functions via Equations S42,S53 and S54.
To parameterize the submodule, density of checkpoint molecules, binding affinity between ligands and receptors, and Hill function parameters were required. 2D densities of CD47 and its receptor on macrophages, SIRPa, were estimated by in vitro assays (Morrissey et al., 2020;Subramanian et al., 2006). Because the absolute number of PD-L1 on cancer cells was not available from literature, PD-L1 density on cancer/ immune cells in the tumor was estimated based on the in vitro measurement on mature dendritic cells and the percentage of tumor/immune cells (45%) that had concurrent PD-1 and PD-L1 expression in TNBC (Cheng et al., 2013;Gatalica et al., 2014;. It was known that CD80 is expressed in TNBC cell lines (Navarrete-Bernal et al., 2020), which could interact with PD-L1 in cis, so the density of CD80 on cancer cells was also estimated by that on mature dendritic cells. Further, binding affinities between all ligands and receptors were available from the literature (Cheng et al., 2013;Jansson et al., 2005). Assuming a Hill coefficient of two for all checkpoint-mediated inhibitions, other parameters, such as the PD-1 density on TAMs and half-saturation constants were estimated to match the experimentally measured phagocytosis activity when one of the checkpoints was absent or blocked by an antibody. Figure 2B shows the overall inhibitory effect of checkpoints based on the checkpoint expression status from model simulation after parameterization. Simulated data points were generated by randomly sampling the density of checkpoint molecules (according to the distributions in Table S1) and calculating the overall inhibitory effect with calibrated Hill function parameters using Equation S54. Negative checkpoint status was simulated by setting the corresponding checkpoint density to 0. Consistent with published experimental studies, the model predicted that phagocytosis of cancer cells by TAMs was about eight times higher when treated by an anti-CD47 antibody (Willingham et al., 2012) and 2-3 times higher when compared with PD-1-negative TAMs (Gordon et al., 2017). Figure 2C shows the dynamics of immune cell subsets in the tumor upon model calibration. Pro-inflammatory macrophages and MDSC are first recruited into the tumor by CCL2, which is followed by the infiltration of CD8 + effector and CD4 + helper T cells and eventually accumulation of regulatory T cells. iScience Article macrophages and CD8 + T cells reach 1e5 cell/mL at day 2 and 14, respectively. Sensitivity analysis shows that the most influential parameters are densities of PD-1 and SIRPa on TAMs, Hill function parameters for inhibition on phagocytosis ( Figure 2D); and rates of TAMs recruitment, polarization, and phagocytosis for tumor growth ( Figure 2E).

Revisiting analysis of the IMpassion130 trial with the macrophage module
After integrating TAMs into the QSP platform, we investigated if the model retains the predictive power on efficacy prediction for clinical trials. To this end, we performed in silico virtual clinical trials using the same approach as our previous analysis of atezolizumab and nab-paclitaxel treatments , but with a recalibrated virtual patient population. We inherited the virtual patient distributions from our previous analysis , which, along with the newly added modules and parameters, were (re)calibrated based on recently published data, such as intratumoral cytokine concentration, density of TAMs, and checkpoint expression. Due to the nonlinear nature of the QSP model, the outputs (e.g., CD8/Treg and CD4/Treg ratios) from simulation with median parameter values did not correspond to the median characteristics of the virtual patient population. Therefore, the descriptive statistics (i.e., median, standard deviation, upper and lower bounds) for parameter distributions were manually adjusted within the physiologically reasonable ranges so that the median characteristics of the virtual patient population matched the clinical measurements. Table 1 lists the model-predicted pre-treatment characteristics of the virtual patient population in comparison with clinically measured values from different sources. 3D cellular densities were either directly reported by clinical studies (e.g., MDSC) or estimated from 2D densities as described above (e.g., TAMs, T cells). Intratumoral cytokine concentrations were measured in a cluster analysis of TNBC biopsy samples (Autenshlyus et al., 2021).
Treatment-related parameters were then recalibrated based on results from a phase 1 clinical trial of atezolizumab monotherapy (Emens et al., 2019) and the placebo comparator arm (placebo + nab-paclitaxel) of the IMpassion130 trial (Schmid et al., 2018). We generated about 1000 virtual patients and conducted in silico virtual clinical trials of atezolizumab and nab-paclitaxel monotherapies using the same dose  Figure 4A shows the Kaplan-Meier curve for model-predicted duration of response on nab-paclitaxel monotherapy. Overall, the model-predicted ORR and DOR are consistent with the clinical results, which fall within the 95% confidence intervals of model predictions. When further characterizing virtual patients into subcategories (i.e., complete/partial response, progressive/stable disease), we observed a shift of more progressive disease to more stable disease in model prediction compared to clinical results. We hypothesized that this  iScience Article shift was because of our assumption on the minimum duration for stable disease (i.e., 8 weeks) and the lack of model prediction of new metastatic lesion, which would be categorized as progressive disease by RECIST criteria (Eisenhauer et al., 2009). As a result, we recategorized the response status in the atezolizumab monotherapy using immune-related response criteria (irRC), by which patients are considered to have partial response or stable disease even if new lesions are present, as long as the respective thresholds of response are met (Wolchok et al., 2009). Figure 3C shows that irRC alleviated the overestimation of stable disease by the model. Interestingly, the model-predicted ORR is lower than the clinically reported value when using irRC, because of its stricter criteria for partial response (i.e., a 50% tumor shrinkage is required), which suggests that more patient data are needed for model calibration to improve its predictive power.
To validate our model structure, we conducted an in silico trial of combination treatment using atezolizumab and nab-paclitaxel with the recalibrated virtual patient population (Table S1). According to the dose regimen in the experimental arm of the IMpassion130 trial, 800 mg atezolizumab was administered every 2 weeks (Q2W) with 100 mg/m 2 nab-paclitaxel administered Q3/4W. Figures 3D and 4Bshow that the model-predicted ORR and DOR for combination treatment of atezolizumab and nab-paclitaxel are consistent with clinically reported values. In our previous analysis, a zero complete response rate was predicted in atezolizumab monotherapy with a limited increase when atezolizumab was combined with nab-paclitaxel . With newly incorporated dynamics of TAMs, the increase in complete response in the combination treatment is more consistent with clinical observation, likely because of the restoration of phagocytosis activity by TAMs upon PD-L1 blockade. To further investigate the role of TAMs in this combination therapy, we performed a subgroup analysis by dividing the virtual patients into subgroups by the medians of model variables of interest and calculating the ORR in each subgroup. Figure 5 suggests that density of TAMs and their subsets, as well as checkpoint expressions on TAMs and cancer cells, may not be predictive biomarkers for this combination regimen, as the confidence intervals overlap between subgroups. As atezolizumab and nab-paclitaxel show additive effect when they are combined in clinical trials and our simulations, majority of the responders in the combination regimen are most likely because of the addition of chemotherapy. As a result, the level of tumor-infiltrating lymphocytes and PD-L1 may not be as predictive as in anti-PD-L1 monotherapy, which relies on restoration of anti-tumoral activity by immune cells through blockade of the immune checkpoint.
We further divided the virtual patients into five subgroups based on the level of each variable of interest and calculated the ORRs in each subgroup. In Figure 6, the differences in ORR were greater than 10% when the subgroups were generated based on density of M2-like macrophages, M1/M2 ratio and expressions of PD-L1, SIRPa, and PD-L2. Specifically, ORR increases as PD-L1 level increases which meets our expectation because a high PD-L1 expression is correlated with high T cell infiltration and thus better Duration of response is defined as the time from the achievement of a response to progression. The median durations of the response with 95% bootstrap confidence intervals are 5.6 (5.6-7.5) and 7.5 (5.6-9.3) months, respectively. iScience Article response to the anti-PD-L1 treatment (Kitano et al., 2017). On the other hand, ORR decreases as SIRPa and PD-L2 densities decrease, suggesting their potential roles in treatment resistance. Surprisingly, ORR decreases as density of M2-like macrophages and M1/M2 ratio increase, even though M2 macrophages exhibit only immunosuppressive activities and have shown a positive correlation with tumor growth without treatment ( Figure 2E). This prediction could potentially result from the correlation between M1-like macrophages and immunosuppressive species, including TGF-b that facilitates polarization to M2 phenotype ( Figure S1). This correlation was also observed in the clinical analysis by (Oshi et al., 2020). Also, the genetic markers used for subtype classification of TAMs differ across studies, which could impact the correlations observed in clinical analyses (Newman et al., 2015;Zhang et al., 2021b). Therefore, there is a need for better understanding of TAMs in the context of treatments (Xiang et al., 2021). Importantly, we have demonstrated here that when virtual patients are calibrated based on population-level data from preclinical studies and clinical trials of monotherapy, the QSP model can make reliable efficacy prediction for the combination treatment.

DISCUSSION
We have here expanded our previously published QSP platform by integrating a macrophage module to investigate the dynamics of TAMs in TNBC. Model calibration, which involves both the new and the preexisting modules as well as the virtual patient distributions, were described in detail above and in STAR methods. In practice, the QSP platform was originally built with seven modules, each of which was written by a MATLAB script with built-in functions of SimBiology Toolbox (Sové et al., 2020). New modules were iScience Article then developed as independent user-defined MATLAB functions that integrate new model elements into the QSP model. In this way, modules can be readily added and removed according to the aims of the study. In fact, the QSP platform has been applied to predict tumor response to various treatments in multiple cancer types with different modules incorporated based on study objectives and data availability for model calibration. This high-usability feature greatly facilitates the wide application of such a platform model in immuno-oncology translational research and drug development (Chelliah et al., 2021).
As discussed in this study, the macrophage module allowed us to investigate the roles of M1-and M2-like TAMs during immunotherapy and chemotherapy in TNBC and improved the ORR predictions. In particular, PD-1 is expressed on both TAMs and effector T cells (Teff) so that both cytotoxic activity of Teff and phagocytosis are inhibited by PD-1-PD-L1/2 interactions, providing an additional target for anti-PD-1/L1 treatments (Gordon et al., 2017). Further, we previously hypothesized that PD-L1 expressed on T cells can interact with CD80 on APC and thus block the co-stimulatory signals (i.e., CD28 À CD80/86) during T cell activation. However, recent findings have suggested that only cis interactions occur between CD80 and PD-L1 (Zhao et al., 2019b). Also, PD-L1 heterodimerizes with CD80 and selectively weakens its interaction with CTLA-4 but not CD28, which provides a mechanistic explanation for the clinically observed synergistic Figure 6. Effects of model variables on response status For each variable, the virtual patient population was sorted by the pre-treatment variable level in ascending order, and evenly divided into five subgroups. The response status of each subgroup in the combination therapy is plotted against the corresponding median variable level. Blue represents partial or complete response. Green represents stable disease. Red represents progressive disease. Related to Figure S1.

OPEN ACCESS
iScience Article effect of anti-PD-L1 and anti-CTLA-4 treatments (Zhao et al., 2019b). In addition to PD-1-mediated inhibition, CD47 also serves as an immune evasion mechanism developed by cancer cells. As shown above, phagocytosis by TAMs can be enhanced eight-fold in breast cancer when CD47 is blocked (Willingham et al., 2012), which suggests that it could be a potential target for cancer treatment, and the present QSP platform may serve as a tool to make prospective efficacy predictions for this emerging strategy and its combination with other checkpoint inhibitors.
In the model, we introduced M1-and M2-like TAMs as two extreme states. In reality, the states constitute a continuum, and we have recently developed systems biology models to describe the complexity of macrophage signaling under different conditions (Zhao and Popel, 2021;Zhao et al., 2019aZhao et al., , 2021. Using these modeling approaches, additional mechanistic detail could be accounted for in extensions of the macrophage module to further enrich the description of macrophage compositions in the tumor microenvironment. It should also be noted that the mechanistic inclusion of TAMs into the QSP platform is one of the several important aspects of the tumor microenvironment, as there are other cell components, such as natural killer (NK) cells, B cells, and fibroblasts that are known to influence cancer progression and treatment efficacy. However, such extensions of the platform have to be carried out step-by-step with careful iterative calibration and validation against cancer-specific experimental and clinical data (Makaryan and Finley, 2020; Makaryan et al., 2020). Although we have made every effort to make modules independent, integrating a module, as demonstrated in this study, is a meticulous process. New model reactions should be detailed enough to reflect the complexity of biological mechanisms that occur in reality, but simplified enough so that the model parameters can be identified by limited data. Once integrated, pre-existing parameters should be recalibrated to make sure that updated model outputs still fall within the physiologically reasonable ranges. This process can lead to a comprehensive mechanistic characterization of the tumor microenvironment that would provide rich therapeutic targets and valuable insights in the basic-to-translational quest in immuno-oncology research.
In silico virtual clinical trials, as we have explored here, is an emerging field of study in immuno-oncology (Chelliah et al., 2021). To capture the inter-individual heterogeneity of patients with metastatic TNBC, we estimated the virtual patient distribution according to pre/clinical data on TNBC or other breast cancer subtypes if TNBC data are not available. In addition, multi-omics and digital pathology data have proven useful when generating physiologically reasonable virtual patients, and in future such data can be further collected from TNBC and incorporated into our QSP platform (Lazarou et al., 2020;Zhang et al., 2021a). Notably, in the tutorial created by (Lazarou et al., 2020), they listed the multi-omics data that can be utilized to predict neoantigen binding affinity with MHC molecules, T cell receptor repertoire diversity, immune cell composition in lymph nodes, and more for QSP model calibration. More importantly, individual tumor growth trajectories from clinical trials could provide substantially more information for virtual patient generation than overall response rates. We also aim to match the simulation setting with the actual clinical trial setting, and we followed the same measurement frequency of tumor size and clinical criteria used in those trials when making efficacy predictions. Nonetheless, the current trial simulation workflow can be further improved by future efforts, as certain factors that lead to treatment discontinuation or categorization of progressive disease, such as appearance of new lesions, severe adverse events, and death, were not described during the in silico trials presented here. For example, unlike in chemotherapy where clinical response often ends by tumor progression led by chemo-resistance, the response in immunotherapy mostly ends when new lesions are found while tumor size is below the response threshold. As a result, the lack of model prediction of new lesions makes the model-predicted duration of response not comparable to the clinically reported value. Additionally, limited tumor detectability of imaging modalities used for tumor size measurement may also affect the response status. To be consistent with our previous analysis, we assume that patients with tumor smaller than 2 mm are categorized as complete responders. These factors should all be considered when interpreting the model predictions, especially for the subcategorization within responders and non-responders.
Although the accelerated approval of atezolizumab for PD-L1-positive metastatic TNBC has been withdrawn since our previous analysis, results from the IMpassion130 study provide useful information for model calibration and validation. The withdrawal of atezolizumab was because of two developments. iScience Article resulting in full regulatory approval for the first-line therapy of patients with metastatic PD-L1-positive TNBC. These results suggest that it may be useful to identify surrogate markers for survival prediction after treatments using the QSP platform, but patient and treatment-specific data are required (Gion et al., 2021). Also, visceral metastasis is correlated with poor survival outcomes in TNBC (Wang et al., 2020b). In the present model, we have only incorporated one tumor compartment that represents an average tumor lesion in metastatic TNBC, with two cancer clones (i.e., sensitive and resistant clones to nab-paclitaxel). Adding multiple tumor compartments that describe the formation and dynamics of metastatic lesions can better assist the study on tumoral heterogeneity and improve efficacy prediction by the model.

Limitations of the study
One of the central factors that could influence the predictive power of in silico clinical trials using a QSP platform is the mechanistic hypotheses that aim to reflect the real biological processes. In our case, hypotheses are made based on our current understanding of the tumor-immune system, which is constantly updated by new experimental data, such as the cis interaction between CD80 and PD-L1 as discussed above. However, some biological processes are still not fully understood, such as the role of PD-L1 on T cells and mechanisms of action (MOA) of PD-1/L1 inhibitors (DeSousa Linhares et al., 2019;Kazanova and Rudd, 2021;Zhao et al., 2019b). This may result in a discrepancy between model prediction and clinical observation. In addition, some simplifications have to be made because of the lack of data during calibration of certain model components. For example, CD4 + helper T cells can be further categorized into multiple subsets (e.g., Th1, Th2, and Th17), each of which may play a unique role in the immune system. However, the lack of quantitative measurements of their cellular density in TNBC prevents them from being explicitly incorporated into the QSP model (Tay et al., 2021). Also, quantitative data were rarely reported for metastatic TNBC lesions, which are critical to study the tumoral heterogeneity (Pasha and Turner, 2021).
Another limitation is related to model parameterization, which has to account for the measurement uncertainty from experimental techniques, cross-study variability, and inter-individual heterogeneity. This is particularly reflected in virtual patient generation. In this study, we manually adjusted the descriptive statistics for each parameter distribution during virtual patient calibration, and randomly generated deviations from the median parameter values that are fitted to the experimental/clinical data. Given that the descriptive statistics and the distribution type can sometimes take more than one fixed set of values, this would have an impact on the model prediction as well as the parameter sensitivity. As a result, larger-than-reality ranges could be generated for parameters and model outputs, and without a sufficiently large dataset of relevant clinical measurements it is hard to apply a robust algorithm to narrow down the virtual patient population while preserving the desirable inter-individual heterogeneity. In addition, covariance between parameters was not considered because of the lack of patient data, as parameters are independently sampled during virtual patient generation, and thus could result in virtual patients that would not have existed in reality (Tivay et al., 2020). In summary, more precise selection criteria and screening methods need to be added into the overall virtual patient generation workflow, as additional data become available, to improve the predictive power of such model-based in silico trials in immuno-oncology research.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

DECLARATION OF INTERESTS
Dr. Emens has had research funding from Genentech, F Hoffman La Roche, EMD Serono, Merck, AstraZeneca, Tempest, Bolt, Silverback, Takeda, CytomX, Compugen, Abbvie, BMS, Next Cure, Immune Onc. She has served as a paid consultant for, F Hoffman La Roche, Genentech, Macrogenics, Lilly, Chug, Silverback, Shionogi, CytomX, GPCR, Immunitas, DNAMx, Gilead, and Mersana. Dr. Emens also has an executive role at the Society for Immunotherapy of Cancer and has ownership interest in Molecuvax. Dr. Popel is a consultant to AsclepiX Therapeutics and CytomX Therapeutics. He receives research funding from AstraZeneca, Boehringer Ingelheim, CytomX Therapeutics. Dr Cesar Santa-Maria has research funding from Pfizer, AstraZeneca, Novartis, Bristol Meyers Squibb and has served on advisory boards for Bristol Meyers Squibb, Genomic Health, Seattle Genetics, Athenex, Halozyme and Polyphor. The terms of these arrangements are being managed by the Johns Hopkins University in accordance with its conflict-of-interest policies. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. Dr. Santa-Maria has research funding from Pfizer, Astrazeneca, Novartis, Bristol Meyers Squibb and has served on advisory boards for Bristol Meyers Squibb, Genomic Health, Seattle Genetics, Athenex, Halozyme and Polyphor.

Lead contact
Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Hanwen Wang.

Materials availability
This study did not generate new materials or reagents.
Data and code availability d All original code has been deposited at Zenodo and is publicly available as of the date of publication. DOIs are listed in the key resources table.
d Any additional information required to reanalyze the data reported in this paper is available from the lead contacton request.

Model overview
The model comprises four compartments (central, peripheral, tumor, and tumor-draining lymph nodes) and ten modules, including a new macrophage module introduced in the present analysis. Each module describes the kinetics and dynamics of a major molecular or cellular species in the tumor microenvironment, such as cancer cells, T cells (i.e., effector T cells, regulatory T cells, helper T cells), and antigenpresenting cells (APCs). All modules from our previous analysis of the IMpassion130 trial (Wang et al., 2021) were inherited, modified, and recalibrated with the new macrophage module, using published data on triple-negative breast cancer (TNBC). The cellular dynamics involved in each module are described below, and the effects of molecular interactions (e.g., immune checkpoints and cytokines) are implemented using

Cancer module
To investigate the dynamics of tumor vasculature and the anti-angiogenic activity of chemotherapy, we adapted the modified Gompertzian model proposed by (Hahnfeldt et al., 1999). Because the tumor volume in our model is calculated by the total number of cancer cells in the tumor compartment, we modified the equations to describe cancer cell dynamics, where the tumor-carrying capacity, K, is translated to the maximal number of cancer cells supported by the tumor vasculature, C max . dC i dt = k C;growth C i log C max C total À k C;death + k C;nabp H nabp min À C total ; K C;nabp Á C total ) dC max dt = k K;g C total c vas c vas + c vas;50 À k K;d C max ðC total V cell Þ 2 3 À k K;nabp C max ½NabP (Equation S2) V T dc vas dt = k vas;Csec C total + k vas;Msec ½Mac M2 À k vas;deg c vas V T + k vas;nabp C total H nabp;vas (Equation S3) Equation S1 describes the dynamics of cancer cells in each cancer clone, C i , added by the cancer module. The first term describes cancer cell proliferation with growth rate constant, k C, growth , total number of cancer cells in all cancer clones, C total , and the maximum number of cancer cells reflecting the tumor carrying capacity, C max . The second term describes the death of cancer cells because of apoptosis, cytotoxic action of nab-paclitaxel, phagocytosis, and effector T cells (Teff). Apoptosis, which is caused by natural cell death, is assumed to be a first-order reaction with the rate constant, k C, death Cytotoxic activity of nab-paclitaxel is incorporated with the killing rate constant, k C, nabp , a Hill function with varying effective concentrations, H nabp , and the number of cancer cells that is accessible by nab-paclitaxel, k C, nabp . Rate of phagocytosis by TAMs depends on the phagocytosis rate constant, k M1, phago , ratio of effector and target cells, and the Hill functions for the inhibitory effects of checkpoint molecules, H Mac, C , and IL10, H IL, phago . Cancer cell killing by Teff is described by the killing rate constant, k C, T , the ratio of Teff and target cells, the ratio of Teff and regulatory T cells ( Equation S2 calculates the maximal number of cancer cells supported by the tumor vasculature. The first term represents the growth of tumor vasculature induced by angiogenic factors, c vas , with a rate constant, k k, g . The second term represents the endogenous inhibition of existing tumor vasculature, such as endothelial cell death, with a rate constant, k k, d (Hahnfeldt et al., 1999). The third term represents the inhibition of tumor vasculature because of nab-paclitaxel, with a dose-dependent rate, k K, nabp (Mollard et al., 2017). The tumor angiogenic factor, c vas , is assumed to be secreted by cancer cell and M2-like macrophages and induced by nab-paclitaxel (Equation S3). The secretion rates are fitted to VEGF-A level measured in preclinical studies (Volk et al., 2008;Wu et al., 2010).
The tumor growth parameters (i.e., k C,growth , k K, g , k K, d , and c vas50 ) and the initial tumor carrying capacity, C max (0), are fitted to the tumor growth curve reported by preclinical studies using TNBC xenograft model. Fitting was performed with the simplex search method using a MATLAB function, fminsearch, to minimize the mean squared difference between observed and predicted values (Lagarias et al., 1998). Using the fitted values as medians, k C, growth and k K, g are varied assuming lognormal and uniform distributions, respectively, to capture the tumor heterogeneity. The median, 60th and 90th percentiles of the simulated tumor growth are plotted with the experimental measurements in Figure S2. k C, growth , k K, g , and k K, d are then scaled up to human using Equation S4. The allometric scaling has been tested in models of various cancer types in human (Elassaiss-Schaap, 2010;Garcia-Cremades et al., 2019;Lindauer et al., 2017;West et al., 2002), while its robustness in translating growth of breast tumor xenografts to human needs to be investigated when clinical data become available.
q human = q mice WT human WT mice iScience Article cancer cell, T cell, and macrophage, respectively. V e is the volume fraction of the intracellular space in breast tumors, which was calculated assuming vascular and interstitial space occupy roughly 2 and 61% of the total tumor volume, respectively (Finley and Popel, 2012).
T cell modules (Teff, Treg, and helper T cell) Naïve T cell dynamics Dynamics of naïve T cells are incorporated into the central (C), peripheral (P), and tumor-draining lymph node (LN) compartments. Equations S6-S8 describe dynamics of naïve CD4 + and CD8 + T cells, where [nT] i represents the average number of naïve CD4 + /CD8 + T cells of a single clonotype in the corresponding compartment.
d dt ½nT C = Q nT ;thym div À Q nT;P;in ½nT C + Q nT ;P;out ½nT P À Q nT ;LN;in ½nT C + Q nT ;LN;out ½nT LN À k nT;death ½nT C (Equation S6) d dt ½nT P = k nT ;pro div ½nT P ½nT P + K nT;pro div + Q nT;P;in ½nT C À Q nT ;P;out ½nT P À k nT;death ½nT P (Equation S7) The initial amount of naïve T cells is calculated by dividing the absolute number of naïve T cells measured from healthy individuals by the T cell clonotype diversity (Autissier et al., 2010;Robins et al., 2009). Q nT, thym represents zero-order thymic export of naïve T cells, whose rate is shown to be correlated with age by (Bains et al., 2009;Ye and Kirschner, 2002) and thus is estimated by the average ages of patients with breast cancer at diagnosis (Yeh et al., 2017). div represents the clonotype diversity, which is 1.11e6 and 1.16e6 for naïve CD8 + and CD4 + T cell, respectively (Robins et al., 2009). Because the naïve T cell densities are sustained mainly by selfrenewal in peripheral lymphoid organs, their proliferation is assumed to occur in the peripheral and the tumordraining lymph node compartments (first terms in S7 and S8) with a rate constant, k nT, pro , estimated based on the in vivo measurements reported by (den Braber et al., 2012). The naïve T cell trafficking among the three compartments is adapted from the model by (Zhu et al., 1996), and the transport rates are estimated to fit naïve T cell densities at the steady state when cancer cells are not present (i.e., no naïve T cell activation by tumor antigens) to the measured naïve T cell levels in healthy individuals (Autissier et al., 2010). When cancer cells are present, Teff is activated from naïve CD8 + T cells, and Treg and helper T cell (Th) are activated from naïve CD4 + T cells. Figure S3 shows that the pre-treatment distribution of naïve T cells (see In silico clinical trial for definition of pre-treatment). The median densities of naïve CD4 + and CD8 + T cells (6.6e5 and 4.7e5 cell/mL) are about 23% and 8% lower than those in healthy individuals (8.6e5 and 5.1e5 cell/mL), which is consistent with the clinical measurements (Autissier et al., 2010;Hueman et al., 2007).

T cell activation and homing
The activation of naïve T cells in the tumor-draining lymph nodes (TDLNs) depends on the number of T cells that can simultaneously interact with mature antigen-presenting cells (mAPCs) and the strength of T cell receptor (TCR)-peptide-MHC (pMHC) interactions, which are implemented as Hill functions named H APC and H Ag , respectively (see the Antigen-presenting cell module and the Antigen module). Equations S9 and S10 describe the dynamics of proliferating T cellsupon activation, iScience Article T cell clones is estimated by the number of neoepitopes identified in TNBC (Morisaki et al., 2021;Narang et al., 2019). k T, act is the maximal activation rate of naïve T cell by mAPCs, k T, pro is the proliferation rate of activated T cell, and N aT is the total number of divisions that an activated T cell can undergo while transitioning to its final form.
In Equation S11, N TCR , N costim , and N IL2 represent the number of cell divisions by signals from TCR, costimuli on mAPCs, and IL-2 secreted by activated CD4 + helper T cells, respectively. According to the experimental data from (Marchingo et al., 2014), the effect of the three signals on activated T cell expansion can be estimated by the linear sum of the underlying signal components.
Similar to naïve T cells in the model, transport of activated T cells (i.e., Teff, Treg, and Th) is adapted from the model by (Zhu et al., 1996) and is described by Equations S12 and S13. The tumor infiltration is limited by the term C 2 total C 2 total + KC;rec to constrain T cell infiltration when the number of cancer cells is low, reflecting the effect of chemoattractant (dePillis et al., 2006). This term is also used to approximate the clearance of immune and myeloid cells because of the lack of pro-inflammatory signals upon tumor eradication in Equations S14, S15, S16, S48, S49 and S50 (Kaech and Cui, 2012;Veglia et al., 2021).
½T C = Q T;LN;out ½T LN À Q T;P;in ½T C + Q T;P;out ½T P À Q T;T;in ½T C C 2 total C 2 total + K C;rec À k T;death ½T C (Equation S12) d dt ½T P = Q T;P;in ½T C À Q T;P;out ½T P À k T;death ½T P (Equation S13) Because of the differential effects of activated T cells on the tumor microenvironment, their equations are specified separately by Equations S14, S15 and S16. For tumor infiltrating Teff, additional death rates are applied to represent Teff inhibition by Treg and cancer cell, which are mediated by IL-10 and PD1-PDL1/2 interactions, respectively (Hsu et al., 2015;Wherry, 2011). For tumor infiltrating Th and Treg, differentiation of Th to Treg is incorporated, which is mediated by TGF-b and arginase-I (Batlle and Massagué , 2019;Serafini et al., 2008).
total C 2 total + K C;rec À k T;death À k Th;Treg H TGFb H ArgI;Treg + k cell;clear K C;rec C 2 total + K C;rec

Antigen-presenting cell module
The APC module describes the APC recruitment into the tumor compartment, APC maturation, and mature APC transport to the tumor-draining lymph node compartment. We assume that the majority of the mAPCs come from the tumor compartment, where they uptake antigens and undergo the maturation process (Lindquist et al., 2004). In Equations S17, S18, S19 and S20, k APC, death is the entry and death rate of APCs; r APC is the baseline APC density; k APC, mat is the maturation rate that depends on concentrations of the maturation signal, IL-12, and inhibitory cytokine, IL-10 (Corinti et al., 2001); k APC, mig is the migration rate of mAPCs from tumor to TDLNs; and km APC, death is the death rate of mAPCs. H APC , which determines the rate of T cell activation in the T cell module, is calculated by Equation S21 based on the type of T cells being activated. That is, n T, clones corresponds to the number of neoantigen clones when calculating H APC for Teff and Th activation from naïve CD8 + and CD4 + T cells, respectively, and it corresponds to the number of self-antigen clones when calculating H APC for Treg activation from naïve CD4 + T cells.

Antigen module
The antigen module is adapted from several well-established models to describe antigen processing and presentation by APC (Agrawal and Linderman, 1996;Chen et al., 2014b;Palsson et al., 2013). Tumor neoantigens and tumor-associated self-antigens are released upon death of cancer cells (Equation S22) and are internalized into intracellular vesicles of APC (Equation S23), where they are processed into short peptides (Equation S24). The peptides then bind with MHC molecules based on their binding affinity (Equation S25), and the pMHC complexes are presented on the cell surface to be recognized by the antigen-specific T cell receptors (Equation S26). Equations S27 and S28 describe the unbound MHC molecules in the endosome and cell surface, respectively. Equations S22, S23, S24, S25, S26, S27 and S28 are applied for both tumor-associated self-antigens and tumor neoantigen clones. P i , p i , M i , and [M P ] i in Equations S23, S24, S25, S26, S27 and S28 represent the average concentration of the antigen, the peptide, the MHC molecule, and the pMHC complexes per APC, respectively. V i and A i represent the volume and the surface area. The subscripts e and s represent the APC endosomal and the surface compartment. The antigen release rate, k dep , is assumed equal to the rate of cancer cell death. k xP, dep is the degradation rate of extracellular antigens. The antigen uptake rate, k up , antigen degradation rate, k P, dep , peptide degradation rate, k p, dep , exocytosis rate of pMHC complexes, k out , and internalization rate of MHC molecules, k in , were estimated by (Chen et al., 2014b). The peptide-MHC association rate, k P, on , was estimated by (Agrawal and Linderman, 1996), and dissociation rate, k P, off , is calculated based on the binding affinity of corresponding antigen, which can vary among different cancer types. For simplification, we assume an average binding affinity for all neoantigen/self-antigen clones to represent the overall effect of all clones.  lick and Renkin, 1970). Using a body surface area of 70 cm 2 /g, the permeability of atezolizumab between the central and the peripheral compartments is calculated to be 2e-8 cm/s, which is used as the starting value for fitting below. Because tumor blood vessels are more permeable than normal capillaries, the permeability between the central and the tumor compartments is estimated to be 3e-7 cm/s according to multiple in vivo studies (Thurber and Wittrup, 2012;Yuan et al., 1995); and the surface area of capillaries per tissue volume is estimated to be 28.4 cm 2 /cm 3 for peripheral tissues and the tumor (Thurber and Wittrup, 2012). As a result, the volumetric flow rates of atezolizumab between the central and the peripheral/tumor compartments are calculated by the corresponding permeability-surface area product. The volumetric flow rate between the central and the tumor-draining lymph node compartment is estimated by an in vivo study (Meijer et al., 2017). a e = 0:483 Ã ðMW Þ 0:386 (Equation S31)

½TCR
Based on the parameter estimation above, we further fit the volumetric flow rate between the central and the peripheral compartments (Q P ), the clearance rate (k cl ), and the volume fractions of plasma (in the central compartment) (g C ) and the interstitial space available to atezolizumab (in the peripheral compartment) (g P ) to match its clinically measured plasma concentration (Stroh et al., 2017). Figure S4 shows the model predicted plasma concentration of atezolizumab with the clinically measured values.
In Equations S32, S33, S34 and S35, [A] i is the antibody concentration, V i is the compartment volume, Q i is the volumetric flow rate between the central and the corresponding compartment, g i is volume fraction of interstitial space available to the antibody, k cl is clearance rate, and Q LD is the rate of lymphatic drainage from tumor to TDLNs and from TDLNs to plasma (Zhu et al., 1996). Subscripts C, P, LN, T represent the central, peripheral, tumor-draining lymph node, and tumor compartments, respectively.

Checkpoint module
Dynamics of PD-1-related checkpoint molecules iScience Article immunological synapse in the present model. We assume that the ligands and receptors are evenly distributed on the cell surface so that their densities in the synapse are calculated by dividing their expression levels by the total cell surface area. We also assume that the explicit representation of the diffusive entry of surface molecules to the synapse is negligible because of its rapid dynamics. Instead, the area of the synapse is increased by a factor of 3 to account for the effect of diffusion (Jansson et al., 2005). The numbers of checkpoint molecules on cell surface are estimated based on measurements using quantitative flow cytometry with fluorescent beads (Cheng et al., 2013), which are then scaled up to account for possible underestimation of PD-L1/2 level by QuantiBRITE bead measurements (Mkrtichyan et al., 2012). These parameters are varied in a wide range in the virtual patient generation to account for the uncertainty and interindividual variability. Equations S36, S37, S38, S39, S40 and S41 describe the PD-1-related dynamics in the model. The model species represent the 2-D densities of the checkpoint molecules in the synapse with a surface area, A syn . k on and k off are the association and dissociation rates of checkpoint interactions; the coefficients for k on and k off of antibody-target binding in Equation S38 represent the stochiometric corrections because of antibody bivalency (Harms et al., 2014); g T is the volume fraction of tumor interstitium that is available to the antibody; c aPDL1 is the intrinsic antibody cross-arm binding efficiency (Harms et al., 2014); the denominator, d syn , is the thickness of the confinement space between two cells during the interaction, which aims to transfer the 3-D binding affinity to 2-D (Jansson et al., 2005); N A is Avogadro's number. k out, PDL1 is the expression rate of PD-L1/2 by IFNg; r PDL1, IFNg is the number of fold increase of PD-L1/2 expression from baseline level by IFNg; k in, PDL1 is the degradation/internalization rate of unbound PD-L1/2 molecules (Shin et al., 2017); r PDL2 is the ratio of PDL2 to PDL1. The number of bound PD-1 molecules to PD-L1/2 molecules in the tumor determines the Hill function, H PD 1, for inhibitory effects of Treg and cancer cell on Teff (Equation S42). iScience Article CD80 and CD86 on APCs. CD28 and CD86 are monovalent, whereas CD80 and CTLA-4 are bivalent. Therefore, CD28 and CD86 can form a monovalent complex, whereas CD28 and CD80, CD86 and CTLA-4 can form bivalent complexes via trans interactions; Also, because of the bivalency of CD80 and CTLA-4, they can form three types of multivalent complexes in a zipper-like fashion (Jansson et al., 2005). In addition, PD-L1 can disrupt CD80 homodimers and form PDL1:CD80 heterodimers via cis interactions, which allows CD80 to bind with CD28 but potentially weaken interactions between CD80 and CTLA-4 (Equations S43, S44, S45 and S46) (Zhao et al., 2019b). CD28 is a co-stimulatory signal that enhances the naïve T cell activation by increasing the number of T cell divisions. Because CTLA-4 outcompetes CD28 because of its higher binding affinity to CD80/CD86, the blockade of CTLA-4 restores ligand availability for CD28 that leads to enhanced T cell activation and proliferation. Similar to those in the PD-1-related dynamics, stoichiometric corrections are incorporated for antibody bivalency and dimerization of CTLA-4 and CD86 on cell surface which also results in bivalency (Bhatia et al., 2005;Harms et al., 2014;Linsley et al., 1995).

(Equation S48)
Overall, the density of checkpoint molecules in the model represents their average expression level on all cancer/immune cells in the tumor, so that the calculated Hill functions represent the overall effect mediated by PD-1 and CD28.

Myeloid-derived suppressor cells (MDSC) module
MDSC module describes MDSC recruitment into the tumor and secretion of arginase-I (Arg-I) and nitric oxide (NO) with their inhibitory effects on T cells. Equation S49 describes MDSC recruitment mediated by CCL2, which is secreted by cancer cells. k MDSC, mig represents the recruitment rate estimated by the median MDSC density in patients with breast cancer (Diaz-Montero et al., 2009). d½MDSC dt = k MDSC;mig V T H CCL2 À k MDSC;death + k cell;clear K C;rec C 2 total + K C;rec ! ½MDSC

(Equation S49)
The major immunosuppressive factors secreted by MDSC are assumed to be Arg-I and NO, whose expression rates are estimated based on in vitro experiments on breast cancer cells . Because only the enzymatic activity of Arg-I is measured in enzyme unit, mU, we use mU as a placeholder of Arg-I concentration in the model, assuming that the protein concentration is proportional to the enzymatic activity. The unit of its production rate is then set to be mU*(microliter)/cell/day to estimate the amount of Arg-I produced by MDSC per day. The units of production rates of NO and CCL2 are set to be ll OPEN ACCESS